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Abstract 

We report on a systematic study of the dynamics of gravitational waves in full 

3D numerical relativity. We find that there exists an interesting regime in the 

parameter space of the wave configurations: a near-linear regime in which the 

amplitude of the wave is low enough that one expects the geometric deviation 

from flat spacetime to be negligible, but nevertheless where nonlinearities can 

excite unstable modes of the Einstein evolution equations causing the metric 

functions to evolve out of control. The implications of this for numerical 

relativity are discussed. 
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Introduction. The study of gravitational waves and their interactions, i.e., the dynam- 
ics of spacetime in its pure vacuum form, is important for both theoretical and observational 
reasons. However, due to its mathematical complexity and the need for previously unavail- 
able large scale computational resources, apart from a few specialized solutions, this area 
of research is largely unexplored. Previous analytic and numerical work on pure gravita- 
tional wave spacetimes, done in ID or 2D, has led to many interesting results These 
discoveries raise interesting questions about waves in 3D. 

Gravitational waves are also about to open up a fundamentally new area of observation: 
gravitational wave astronomy. A new generation of detectors [|| should see waves for the first 
time near the turn of the century. Even though the observed waves are expected to be weak 
here on Earth, they will have been generated in regions with strong, highly dynamic and 
nonlinear gravitational fields. It is therefore essential to study accurately both the strong 
and weak field regimes, as well as the long term secular behavior in transitory intermediate 
regimes. 

Thus we have been motivated to undertake a systematic study of gravitational waves in 
3D numerical relativity over the past few years. 3D studies of general gravitational wave 
spacetimes are different from lower dimensional studies in two very important respects. 
First, unlike in lower dimensional studies, there are no assumed symmetries to fix the gauge. 
One has to consider the full set of gauge freedoms during the numerical evolution. Second, 
3D calculations are bound to have lower resolution, and hence larger errors. This in turn 
makes it difficult to separate physical from numerical and gauge effects. These two effects 
are in fact related, as finite differencing errors tend to excite all gauge freedoms available. 
The excitation of unstable modes, whether they correspond to gauge freedoms or not, is the 
underlying principle of the phenomena discussed in this paper. 

We have developed advanced 3D Cartesian numerical codes, using different numerical 
techniques as well as different formulations of the Einstein equations, to study gravitational 
waves. Full details of our methods and codes will be reported elsewhere In this paper 
we focus on one important lesson we learned with the large amount of numerical studies 
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carried out using these codes. 

When the amphtude of the gravitational waves is low, the evolution is linear. At high 
amplitudes there are nonlinear geometric effects: the energy-momentum carried by the waves 
themselves changes the geometry of the spacetime so that the evolution cannot be described 
by linear equations. We have seen both of these regimes in our studies, as expected. What 
was not anticipated is that beyond the linear regime, but before the amplitude is large enough 
for nonlinear geometric effects to be significant, there exists a regime where the evolution of 
waves, as described by the full nonlinear Einstein equations, can differ dramatically from 
the evolution described by the linearized versions of the Einstein equations. This is what we 
called the "Near-Linear Regime" (NLR). The aim of this paper is to draw attention to the 
existence of such a regime. 

Many questions immediately come to mind: In what way do these effects show up, how 
do they arise, and can they be controlled? What are the boundaries between the linear, 
near-linear, and nonlinear regimes? What are the implications in the numerical study of 
gravitational waves, and gravitational systems in general? Will such effects be important 
for the numerical simulation of waves as weak as those we expect to observe astrophysically? 
These are the questions we address in this paper. 

Statement of the problem. The evolution of the metric functions is described by 
the Einstein equations in the standard 3+1 form 

dtgij = -2aKij + + Vj(3i 

dtK,, = -V.Vja + a [r,, + K K,, - 2K,^K^) (1) 

The question is: When will this evolution be different from what is described by the linearized 
evolution equations, namely, by taking gij = 6ij + ehij, and keeping only terms to first 
order in e in Eqs.(|l|)? More precisely, we ask when the discrete versions of Eqs. (|I]) yield 
results significantly different from the discrete versions of the linearized equations, using the 
same finite differencing techniques. We use the same initial data, which satisfies the initial 
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constraint equations, for both the nonhnear and hnear evolutions. 

The purpose of this paper is to demonstrate, with concrete examples, that significant 
nonlinear effects appear before one might expect. Without careful analysis, these effects 
can be mistaken for true geometric effects, e.g, the formation of a spacetime singularity. 
These effects are due to unstable nonlinear modes admitted by Eqs. (|I]) (and/or their finite 
differenced versions), excitable by waves even with quite low amplitudes. 

Quadrupole waves. The first example is of quadrupole waves similar to those studied 
in [^, except that we use an imploding-exploding wave configuration which is time symmetric 
0, and we take the metric given in |Q as the conformal part of the metric function. We 
then determine the conformal factor using the York decomposition to solve the hamiltonian 
constraint equation 0. The time symmetry of this imploding-exploding data automatically 
guarantees that the momentum constraint is also satisfied. 

We study a case where the wave is initially centered at the origin, with a peak amplitude 
of y4 = Qij — Sij ~ 10~^, and a width of order unity (we use c = G = 1 throughout this paper). 
The evolution of this initial data is predictable: As the amplitude of the wave is of the order 
A ~ 10-3, and both space and time derivatives can bring in only factors of order unity, the 
energy density (e.g. Landau-Lifshitz pseudo energy density) must be of order ~ 10"^. 
With this "wave energy density" concentrated in a spatial dimension of order unity, the 
strength of the "gravitational potential M/R" (or spacetime curvature), can at most be of 
order 10"^. This is approximately 10^ times too weak to produce a significant effect on the 
propagation of the wave, e.g. to "bend" the null characteristics of the spacetime. Hence 
we expect the linear wave to expand outward just as in flat space. The scattering of the 
wave by itself, and any tail effects would affect the metric functions to at most of order 

~ 10^^, which is beyond the level of accuracy in our simulations. Hence the wave train 
should propagate outward basically in a linear fashion, except for a slight delay in phase. 

This picture is confirmed in Fig. |I^ where we show the evolution of a typical metric 
function g^^ (solid lines) along the z axis from the initial time t = to t = 3, with zero 
shift and unit lapse. This is to be compared to what is obtained with the same initial data, 
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but evolved with the hnearized evolution equations (dotted lines). At early times, the two 
waveforms are completely indistinguishable. 

Next we look at the evolution at late times in Fig. |I|b and zoom into the origin where we 
expect the spacetime to have returned to fiat space. The nonlinear results (solid lines) are 
basically one everywhere except for a dip near the origin. This is clearly due to nonlinear 
effects, as there is absolutely no dipping in the linearized treatment (dotted lines). There 
are similar dips in Qyy and Qzz, while the off diagonal metric functions develop a shear-like 
structure. These "residual" features will keep growing in time until the code crashes. This 
behavior is independent of grid resolution and has a characteristic growth rate that increases 
in proportion to A^. 

What causes this result? There are four possibilities: (i) nonlinear spacetime geometric 
effects which are not captured by the order of magnitude estimate given above; (ii) nonlin- 
earity coupled with finite differencing errors; (iii) nonlinearity coupled with gauge effects; 
(iv) bugs in the nonlinear evolution code. In principle, any of these possibilities can cause 
this feature. 

After much investigation that includes numerous code tests, convergence studies, 
comparisons with linear evolutions, different numerical techniques, and even different for- 
mulations of the Einstein equations, we conclude that, in this case, it is (iii). The nonlinearity 
excites an unstable gauge mode of the evolution equations. A posteriori, the effect is simple 
to understand. The potential well created by the initial wave packet, although quite small 
in amplitude, sets the coordinate lines into motion, drifting towards each other. As the wave 
propagates outward, the potential well (spacetime curvature) reduces to zero, but the coor- 
dinate lines, once set into motion, will keep drifting towards the center. Hence over time, 
the distance between, say, constant x coordinate lines decreases and a dip in g^x develops. 

At what point does this effect appear (the lower boundary of the NLR)? This depends 
on the amplitude of the wave and the length of time it is evolved. We have observed dipping 
for waves with magnitude A ~ 10^^ and spatial extension the same as in Fig. |I|a, after it is 
evolved to about t = 10. 
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At what amplitude would nonlinear geometric effects become significant (the upper 
boundary of the NLR) for the quadrupole waves at late time? In Fig. |l|c we show the 
Riemann invariant J ^ obtained using nonlinear and linear evolutions at time t = 4, for 
various initial wave amplitudes A's. We see that despite the dipping of the metric, the 
invariants are basically the same up to A ~ 0.05. Beyond that, the linear and nonlinear 
results rapidly diverge, indicating nonlinear geometric effects. 

The existence of unstable modes in the nonlinear evolution equations in the NLR high- 
lights the importance of controlling the motion of the coordinates. In Fig. |l]d we show the 
evolution of Qxx at various times using the minimal distortion shift and unit lapse. We 
see that there is no dipping. Details of this calculation and other methods in controlling the 
coordinate motion with both shift and lapse will be discussed elsewhere 0]. 

Colliding Wave Packets. Next we discuss the collision of two plane wave packets. In 
this case a completely different type of unstable mode is excited in the numerical evolution 
of Eqs. ([^). The initial data is given by two nearly gaussian packets of the form shown by 
the solid line at t = in Fig. ^ for the metric function g^x- The evolution is carried out 
with a = 1, zero shift, and grid spacing Ax = 0.05. The wave packets, which are initially 
located aX z = ±3, travel towards each other and collide at the center of the grid at t = 3. In 
the wake of the waves, following the collision, we see that the metric function Qxx develops 
an upward drift. Comparing gxx from the linear (dashed line) and nonlinear evolutions at 
t = 9, we see clearly that the drift is a nonlinear effect. 

Again we face the four possibilities pointed out above. In fact it is known that when 
two plane symmetric waves collide, a curvature singularity is generated due to the focusing 
effects of the waves even for arbitrarily weak waves. However, it is easy to show, based on 
the colliding packet study in 0, that the singularity will develop at a time t ~ A^/ [(271^4) ^a] 
after the collision, where A is the characteristic wavelength, a is the characteristic width of 
the packet, and A is the characteristic amplitude of the packet. For the case here, with 
A ~ 1, 0" ~ 1, and A ~ 10~^, we expect the singularity to appear at t ~ 250, which is far 
beyond any evolutions shown here. Of course this estimate does not rule out that this drift 
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is a "precursor" of the singularity, nor the possibihty of other geometric nonhnear effects. 

After much investigation we confirm that the drift is an unstable mode in the nonlinear 
evolution equations. (Q). We note that this is not in contradiction to the expectation that 
the Einstein equations are stable for weak waves (weak perturbations of the flat spacetime). 
It is the constraint equations that rule out these unstable modes. In our free evolution code 
the constraint equations are not enforced. This allows the unstable modes to develop after 
they are excited by the numerical errors in the evolution. The amount it is excited depends 
on the details of the numerical scheme. The instability causes the metric components in 
the wake of the waves to drift upwards according to gij — 6ij oc — ln(Ct + 1). The time 
scale 1/C of the instability depends strongly on the finite differencing scheme used. For the 
same run parameters used in Fig. 0, we find that C = — 1.6 x 10"'' for a staggered leapfrog 
scheme, and a very different C = 7.9 x 10^^ for a MacCormack scheme. We stress that the 
unstable mode exists on the level of the evolution equation and is excited in both of these 
two standard schemes, although with very different rates. 

Conclusion. We have shown that there exists a "near-linear regime" (NLR) in the 
numerical study of waves in general relativity, in which the amplitude of the wave is weak 
enough so that nonlinear geometric effects are unimportant, but nevertheless unstable modes 
of the Einstein evolution equations can be excited by the nonlinearity. The instabilities 
manifest themselves as secular changes in the metric functions, which eventually grow to a 
level that affects the evolution, even crashing the code. There are various types of unstable 



modes. We have shown two examples of such modes [|10|. In the quadrupole wave case, 
the instability is caused by nonlinearity coupled with gauge freedom. In the colliding plane 
wave packets case, the instability is caused by nonlinearity coupled with finite-differencing 
errors. The former effect conserves the constraint equations, the latter does not. In general, 
effects will be cross coupled and show up in many different forms. Further details of this 
analysis, as well as studies of other gravitational wave spacetimes will be given elsewhere. 
Our studies have been based solely on free evolutions. Unstable modes will be different for 
constrained evolutions. 



The aim of this paper is to draw attention to the existence of the NLR, so that the 
large amount of work leading to our keen awareness of these NLR phenomena need not 
be duplicated by other research groups in the field. We believe the instabilities in this 
regime are important in general for numerical relativity, with implications not just for grav- 
itational waves, but also for evolving other 3D gravitational systems. The determining 
factor of whether nonlinear effects are important or not depends on both the strength of the 
gravitational potential and on the length of time the general relativistic system is evolved. 
To determine gravitational waveforms for astrophysical events, e.g., inspiral coalescences 
of compact binaries, which are expected to be one of the important observable sources of 
gravitational radiation, we would like to be able to evolve the system, compact objects and 
waves, to hundreds or even thousands of M, where M is the mass of the system. For the 
success of such long time scale simulations, a full understanding and control of the kind of 
unstable modes pointed out in this paper are indispensable. 

This research is supported by NCSA, PSC, and NSF grants Nos. PHY94-04788, PHY94- 
07882, ASC95-03978 and ASC93-18152. 
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FIGURES 

FIG. 1. (a) The metric function Qxx is plotted along the z axis for a quadrupole wave. The wave 
peaks at the origin at i = 0, and disperses outward. By t = 3 the region near the origin returns 
basically to flat space, as expected. The same initial data is evolved with the nonlinear evolution 
equations (solid lines), and the linear evolution equations (dotted lines), (b) Qxx is shown at later 
times t = 4,6,8, 10 for both linear (dotted lines) and nonlinear evolutions (solid lines). Note the 
development of a dip near the origin with the nonlinear evolution, (c) The Riemann invariant J 
at the origin, is plotted against the initial amplitude of the wave Qzz — 1, for both the nonlinear 
(solid line) and linear evolutions (dotted line) at t = 4. At small amplitudes the two J's are the 
same, despite the dipping of the metric in the non-linear evolution. For perturbations of order 5% 
we start to see deviations, indicating the presence of nonlinear geometric effects, (d) An evolution 
of the same quadrupole data, but now evolved with the minimal distortion shift and a = 1. Note 
that the shift is able to hold the coordinate lines from drifting in. 

FIG. 2. The metric function gxx is shown for colliding plane wave packets. The packets are 
initially centered at z = ±3 at t = and are moving towards each other. They collide at the origin 
at t = 3, then pass through each other. The solid lines show the nonlinear evolution. Note the 
upward drift in the "wake" regions at t = 9. Such drifting is absent in the linear evolution (dashed 
line) . 
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